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I.  INTRODUCTION 


Studies  in  separating  boundary  layers  leading  to  accurate  methods  for 
predicting  flow  separation,  predicting  flows  on  the  verge  of  separation, 
predicting  separating  and  reattaching  flows  and  predicting  fully-separated 
flows  find  a large  number  of  applications  in  aerodynamic  problems  and  naturally 
in  the  design  of  both  commercial  and  military  airplanes,  for  example,  an 
accurate  determination  of  the  separation  point  is  very  crucial  in  many  problems 
since  separation  strongly  influences  the  performance  of  acrodv nami c configu- 
rations. A common  procedure  employed  to  determine  this  point  is  to  <^olve  the 
governing  boundary- layer  equations  for  a given  external  velocity  distribution 
and  find  the  point  (if  any)  where  the  wall  shear  goes  to  zero.  According  to 
recent  studies  (for  example,  ref.  i),  such  a prediction  can  be  done  accurately 
for  two-dimensional  and  axi symmetri c laminar  and  turbulent  flows.  Such  a 
capability  should  also  exist  shot cly  for  three-dimensional  flows  in  view  of  the 
considerable  work  being  done  in  this  area  (see  for  example,  refs.  2,  3,  4). 

Predicting  flows  on  the  verge  of  separation  (that  is,  flows  with  ze*-o  wall 
shear)  is  quite  important  in  many  problems.  In  the  case  of  airfoils  where  it 
is  desired  to  maximize  the  lift,  it  is  necessary  to  compute  tho  minimum  distance 
over  which  a given  pressure  rise  can  be  obtained  without  the  flow  separating. 

The  most  rapid  pressure  rise  that  it  is  possible  to  obtain  occurs  when  the  wall 
shear  stress  along  the  suction  side  of  the  airfoil  decreases  to  zero.  There- 
fore, it  is  of  considerable  interest  to  be  able  to  calculate  boundary  layers 
with  specified  values  for  the  wall  shear  that  aUo  decrease  to  zero.  The 
Liebeck  airfoils  discussed  in  refs.  5 and  6 are  designed  on  that  principle. 

Prediction  of  partially  separating  flows  is  one  of  the  most  difficult, 
yet  rewarding  tasks  in  aerodynamics.  An  economic  and  efficient  operation  of 
aerodynamic  devices  depends  on  smooth,  streamlined  flow.  The  upper  limit  of 
this  efficient  operating  range  is  marked  by  the  flow  breakaway,  called  sepa- 
ration, or  stall.  The  ideal  attached-flow  conditions  are  seldom  attained  in 
practice  since  a design  is  a series  of  compromises  between  conflicting  require- 
ments. As  a specific  example,  consider  an  airplane  wing  whose  maximum  lift  is 
determined  by  separation.  If  thi:.  wing  was  designed  from  the  performance 
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standpoint  alone,  simultaneous  spar, wise  stall  would  be  desirable.  However, 
because  nobody  wants  to  stall  without  lateral  control  (by  ailerons),  wings  are 
designed  to  have  progressive  stall  from  the  wing  root  out.  Operating  with 
partially  separated  wings,  especially  for  swept  wings,  creates  in  turn, 
longitudinal  control  problems  because  of  shifts  in  the  center  of  pressure. 

From  the  standpoint  of  design  of  control  surfaces  it  is  necessary  to  have 
methods  for  calculating  overall  forces  on  wings  with  partial  separation.  Today 
such  a capability  for  calculating  partially  separating  flows,  even  in  two 
dimensions,  does  not  exist. 

Prediction  of  separating  and  reattaching  flows  is  also  very  important. 

Two  typical  examples  of  such  flows  are  a leading  edge  bubble  and  shock  boundary- 
layer  interaction.  In  both  cases  the  boundary  layer  separates  from  the  surface 
but  reattaches  after  a short  distance.  Here  it  is  important  to  know  under 
what  conditions  the  boundary  layer  reattaches  or  separates  completely.  The 
reattachment  after  a leading  edge  bubble  on  highly  swept  wings  is  relevant  to 
the  leading-edge  vortex  formation,  and  the  absence  of  reattachment  after  a 
shock  signifies  shock-induced  stall.  Today  a satisfactory  prediction  of 
separating  and  reattaching  flows,  even  for  two  dimensions,  again  does  not 
exist. 

In  recent  years  a number  of  studies  on  inverse  boundary-layer  flows  have 
been  conducted.  Except  for  the  integral -method  approaches,  most  of  these  studies 
have  been  directed  to  laminar  layers.  In  ref.  7,  Catherall  and  Mangier  solved 
the  laminar  boundary-layer  equations  in  the  usual  way  until  the  separation  point 
was  approached.  By  assuming  that  the  displacement  thickness  behaves  in  a regular 
prescribed  fashion  in  the  region  of  the  separation  point,  they  calculated  the 
pressure  distribution  in  that  region  for  the  prescribed  displacement  thickness 
distribution.  Their  numerical  solutions  did  not  show  any  signs  of  a singular 
behavior  at  separation. 

In  ref.  B,  Keller  and  Cebeci  solved  the  laminar  boundary-layer  equations 
for  a prescribed  positive  wall  shear  and  in  ref.  9,  Klineberg  and  Steger  solved 
them  for  a prescribed  negative  wall  shear.  In  ref.  10,  Carter  presented  numer- 
ical solutions  of  the  laminar  boundary-layer  equations  involving  separation  and 
reattachment.  He  obtained  solutions  with  an  inverse  procedure  in  which  he 
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prescribed  the  displacement  thickness  or  the  wall  shear.  He  compared  his 
results  with  Klineberg  and  Steger's  separated  boundary-layer  calculations^ 
and  with  Briley's  solution^^  of  the  Navier-Stokes  equations  for  a separated 
region. 

In  ref.  12,  Cebaci , Berkant,  Si  1 i vri  and  Keller  solved  the  turbulent 
boundary- layer  equations  for  a prescribed  positive  wall  shear.  The  only  other 
turbulent  boundary- layer  calculations  for  flows  with  prescribed  wall  shear  were 
made  by  Kuhn  and  Nielsen'  , by  using  an  integral  technique.  However,  unlike 
ref.  12,  their  solutions  include  negative  wall  shear  as  well  as  positive  wall 
shear. 

The  work  described  in  this  report  is  one  phase  of  studies  on  separating 
flows  conducted  under  the  contract  N00014-74-A-0203-0001 , NR215-234,  from  the 
Of;  ice  of  Nav^l  Research.  It  deals  with  the  calculation  of  laminar  and 
turbulent  boundary-layer  flows  for  standard  and  inverse  boundary-value  problems, 
and  is  applicable  to  both  incompressible  and  compressible  flows.  The  standard 
boundary-layer  problem  considers  the  solution  of  the  usual  boundary- layer 
equations  for  a given  external  velocity  distribution.  The  inverse  problem 
considers  the  solution  of  the  governing  equations  for  assigned  wall  shear  or 
for  assigned  displacement  thickness.  It  provides  a very  useful  and  powerful 
method  for  calculating  flows  on  the  verge  of  separating. 

The  method,  which  is  developed  for  two-dimensional  flows,  can  easily  be 
extended  to  axisymmetric  flows.  It  also  has  the  potential  to  be  used  in  a 
number  of  problems  that  require  inverse  boundary-layer  procedures.  Some  of 
them  are: 

1.  Laminar  flow  control  studies.  Here^the  problem  is  to  find  the  minimum 
suction  rate  to  keep  the  flow  laminar. 

2.  Design  of  ducts  for  a given  pressure  distribution. 

3.  Design  of  optimum  ducts. 

4.  Calculation  of  attached  duct-flows  (inviscid  and  viscous)  such  as 
those  in  diffusers. 

5.  Design  of  two-dimensional  and  axisymmetric  shapes. 

6.  Possible  application  to  separated  external  flows. 

7.  Possible  application  to  separated  flows  in  ducts,  i.e.,  diffusers. 
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II.  GOVERNING  EQUATIONS 


2.1  Boundary-Layer  Equations 

The  governing  boundary- layer  equations  for  steady,  two-dimensional, 
compressible,  laminar  and  turbulent  boundary  layers  are  the  continuity, 
momentum  and  energy  equations.  They  are  given  by: 

Continuity 


lx  -° 


(2.1.1) 


Momentum 


au  . — au 
pu  a7+pvay 


+ fy  [u  ly  — PUrVr]  (2-1.2) 


Energy 


„ 9H  x T 9H  - 3 fy  3H  — rrrr  , /,  1 \ 3ul 

pu  3x  pV  ay  ay  [pr  ay  pv  H y\1  Pr)  u ayj 

Here  -puTvr  and  -pv'H'  denote  the  Reynolds  stresses  and 


(2.1.3) 


pv  = pV  + p'v1 

The  boundary  conditions  for  (2.1.1)  to  (2.1.3)  for  zero  mass  transfer  are: 

y = 0 u,v  = 0 H = Hw  or  (aH/ay)w  = given 

(2.1.4a) 


y - 


u -*•  ue(x)  H ->  He 


(2.1.4b) 


2.2  Closure  Assumptions  for  the  Reynolds  Stresses 

The  solution  of  the  system  given  by  (2.1.1)  to  (2.1.4)  requires  closure 
assumptions  for  the  Reynolds  stresses,  -pu1 v'  and  -pv'H'.  In  our  study 
we  use  eddy  viscosity  (em)  and  eddv  conductivity  (e^)  concepts  and  define 


~pU ' V1  = PE. 


au 
m ay 


. . m i _ aH 
■pVTT  - 0Eh  — 


(2.2.1) 
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and  relate  t and  t,,  to  a turbulent  Prandtl  number  Pr.  by 
m H t 

p-'t  * V‘-H 

/ 1 A \ 

According  to  the  eddy-viscosity  formulation  of  Cebeci  and  Smith'  , tii  ? 
turbulent  boundary  layer  is  divided  into  two  regions,  called  inner  and  outer 
regions,  and  the  eddy  viscosity  is  defined  by  separate  formulas  in  each 
region.  They  are: 


£ 


m 


<s>i  ■ L 


a u | 
ay  | 


(£m}i  - ( Lm}o 


(2.2.3a) 


/ (ue  -u)dy|  . » ('Jo 

0 


(2.2.3bA 


Here 


L - ry  [1  - exp(-y/A)] 


N = [i  -H.8(uw/ue)(oe/pw)2  p+]1/2 


(2.2.4) 


+ 

P 


T 


u 


1/2 


In  (2.2.3b)  and  (2.2.4)  a,  k and  A+  are  "universal"  constants  equal  to 
0.0168,  0.40,  and  26,  respectively,  for  hign  Reynolds  number  flows,  Rfi  > 5000. 

To  compute  flows  at  low  Reynolds  numbers,  one  can  modify  these  three  coefficients 
as  discussed  in  ref.  3,  p.  221. 


The  eddy-- vi scosi ty  formulas  (2.2.3)  can  also  be  modified  to  compute 
transitional  boundary  layers  as  well  as  boundary  layers  in  which  the  stream- 
wise  wall  curvature  becomes  important.  Again  for  a detailed  discussion  see 
ref.  3,  p.  232. 
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2 . 3 Transformation  of  the  Governing  Equations 

Before  we  solve  the  system  given  by  (2.1.1)  through  (2.1.4)  with  the 
Reynolds  stresses  replaced  by  (2.2.1),  we  introduce  the  Falkner-Skan  transforma- 
tion to  remove  the  singularity  at  x = u anu  to  stretch  the  coordinate 
in  the  y direction;  we  define: 


x = x 


(2.3.1) 


This  transformation  not  only  alio  's  the  calculations  to  be  started  very  easily 
at  the  leading  edge  or  at  the  stagnation  point  but  also  removes  the  large  vari- 
ation in  boundary-layer  thickness  along  the  surface.  In  transformed  variables, 
the  velocity  profiles  and  temperature  profiles  do  not  change  "much"  as  the 
calculations  proceed  in  the  x-direction.  'his  results  in  small  computation 
times,  and  allows  larger  spa^ings  to  be  taken  in  the  x-direction. 


We  nexr  define  the  stream  function  g.  by 


D ip 

pu  = *7  • 


PV  = - 


9 

dx 


(2.3.2) 


and  a dimensionless  stream  function  f(x,n)  by 


<*’  = (peueuex)1/,2f(x,n) 


(2.3.3) 


With  the  relations  (2.3.1)  through  (2.3.3),  we  can  write  the  momentum  and  energy 
equations  as 


Momentum 


(bf")'  + P^f"  4 P2[c  - (f)2]  = x (f1  §£--  f" 

Energy 

W * dff-r  * P,fg'  -X  (f  |a-9'  ft 
Here  primes  denote  differentiation  with  respect  to  n and 


(2.3.4) 


(2.3.5) 


f'  =-  , 

ue 


g = h-. 


C s •JE^- 
pe‘e 


(2.3. 6a 1 
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Momentum 
n - 

n = <*> 

Energy 

n = 0 

r,  = n 


0C-  e ‘ Fr  ( 

1 + + fr_\ 

fm  Prt  / ’ 

u2 , 

~^r)  (2.3.6b) 

+ P2  * v7  3*  (pe"e)] 

P~  = — 
2 ue 

du 
e 

dx 

(2.3.6c) 

boundary  conditions  in 

(2.1.4)  become 

0 f ! 

= 0 

f1  - 0 

(2.3.7a) 

f'  = 1 


g = g 
g = i 


w 


or 


= given 


(2.3.7b) 

(2.3.8a) 

(2.3.8b) 


In  terms  of  transformed  variables,  the  inner  and  outer  eddy  viscosity  formulas 
in  (2.2.3)  can  be  written  in  dimensionless  form  as 


<4>i  (~f  “I"  V [t  - exp(-y/A)]‘ 

l/> 


(1  - f )d-:| 


(2.3.9a) 


(2.3.9b) 


where  em  = e„/v  and 
m m 


^J/2|UoWduJ/2  j1/4^„)1/2n/a+ 


t- <''*{?) (rB  ^ 


**¥*•  ''jr'i  (2-3J0) 


u f 

Iw  w 

TTo 

ue  R 


e '^w 
„ \l/2 


vf  u. 

— 7 

uc 


R = 

X v 


U X 

e 
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III.  STANDARD  AND  INVERSE  PROBLEMS 


The  system  given  by  (2,3.4),  (2.3.5),  (2.3.7)  and  (2.3.8)  with  specified 
ug(x)  or  P2 ( x ) is  the  typical  two-dimensional  boundary-layer  problem  for 
laminar  and  turbulent  flows.  For  convenience  we  shall  call  it  the  standard 
problem. 


There  are  a number  of  problems  that  require  inverse  procedures  in  viscous 
flows.  One  type  of  inverse  problem  results  from  requiring  that  the  local  sk in- 
friction coefficient  defined  by 


c 


f 


0/2)ne< 


(3.1.1) 


be  specified.  In  terms  of  transformed  variables  (3.1.1)  becomes 


2f"C 
w w 


(3.1.2) 


Another  type  of  inverse  problem  results  from  requiring  that  the  displace- 
ment thickness  defined  by 


t 

f 


00 


be  specified.  In  terms  of  transformed  variables  (3.1.3)  becomes 


6* 


(3.1.3) 


(3.1.4) 


Other  inverse  problems  can  be  formulated  for  the  problems  discussed  at  the  end 
of  Section  I.  The  system  given  by  (2.3.4),  (2.3.5),  (2.3.7),  (2.3.8)  and  (3.1.2) 
or  (3.1.4)  is  overdetermined  and  we  cannot  specify  P^x)  (i.e.,  up(x)) 
arbitrarily.  Rather  we  must  determine  (x ) as  well  as  f(x,n)  to  solve 

the  system.  In  our  study  we  use  Newton's  method  and  determine  the  unknowns  by 
the  procedure  discussed  in  the  next  section. 
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3.1  Newton's  Method  for  the  Inverse  Problem 


To  describe  our  numerical  approach  to  the  problem  of  the  speci fi ed 
cf  case,  let  us  assume  that  at  x = xn_j  we  are  given  the  profiles  of 
77  f‘,  f",  g,  g',  the  pressure  gradient  P2^xn-1^  ancl  the  ve^oc1ty  ue^xn-l^’ 
At  x = we  seek  an  accurate  approximation  to  the  solution  of  (2.3.4), 
(2.3.5)  subject  to  (2.3.7),  (2.3.8)  for  a given  cf(x).  To  start  the  calcula- 
tions, it  is  necessary  to  know  ?2(x)  and  ug(x).  The  latter  is  necessary 
since  Is  a function  of  ug.  In  our  method  we  assume  P2 ( x ) and  calculate 
u (x)  from  the  definition  of  P-(x)  in  (2.3.6c).  Using  central  differences 
we  approximate  P2  and  solve  it  for  u^  = ug(xn)  to  get 


ph-1/2  7 

n n-1  P2  + 2fln 
ue  ■ 'ue  


- 2c 


n 


where 


(3.1.5) 


n-1/2 


n-1 


P"']/2  = \ (p2  + P^1). 


n-1/2 


2 (xn  + xn-l ^ 


(3.1.6) 


Once  P2(x)  and  u&(x)  are  known,  then  the  standard  problem  (2.3.4),  (2.3.5) 
subject  to  (2.3.7),  (2.3.8)  can  be  solved.  The  numerical  method  used  to  do  this 
will  be  described  in  Section  IV.  Let  us  denote  the  solution  of  the  standard 
problem  by 

(3.1.7) 


f(x,n)  = B[x,n,P2(x)] 


Using  this  solution,  we  can  now  calculate  c^  (which  we  shall  denote  by  Cfc) 
from  (3.1.2).  Recalling  that  the  desired  value  for  the  skin-friction  coefficient 
i s cf(x) , we  form: 


9[Po(x)]  = cf  - cf  (3.1.8) 

c 

and  seek  P2(x)  such  that  >>^^2 )^ = ^ on  x > 0. 


To  solve  <t[P2(*)3  = 0,  we  use  Newton's  method.  With  some  estimate 
P^o)(x)  of  the  desi-ed  pressure  gradient,  we  define  the  sequence  P^(x)  by 
setti ng 
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(3.1.9) 


P^+1(x)  = Pj(x) 


$[P^}(x)] 

0/3P2){^[P^v)(x)]} 


The  derivative  of  <p  with  respect  to  P^  can  be  obtained  from  (?.1.8)  by 
making  use  of  the  relation  given  by  (3.1.2)  and  (3,‘'.r).  This  gives 


where 


2C  1 

u"-1 

2a 

) 

w 

C‘l  4.  ft  1 

n 

(P2"1/2  - 2an)2 

1 

F" 


w 


(3.1.10) 


(3.1.11) 


To  summarize  one  step  of  iteration  of  Newton's  method,  we  first  estimate  a 
value  for  P2(xn),  then  calculate  from  (3.1.5),  and  obtain  a solution  of 
(2.3.4),  (2.3.5)  subject  to  (2.3.7)  and  (2.3.8).  The  solution  yields  a Cf^ 
according  to  (3.1.2).  From  this  result  and  from  the  desired  value  cf(xn),  we 
find  4.  from  (3.1.8).  It  is  then  clear  that  the  next  value  of  P2(xn) 
can  be  calculated  from  (3.1.9),  provided  that  3^/3 is  known.  In  Section  4.3 
we  shall  discuss  its  calculation.  The  iteration  process  is  repeated  until 


|p(v+l  ) 

I 2 


(3.1.12) 


where  y-j  is  a small  error  tolerance. 


Our  procedure  for  the  specified  6*  case,  is  similar  to  the  procedure  for 
the  specified  cf  case.  The  difference  in  the  procedure  starts  after  we  get 
the  solution  ^3.1.7).  uehutin^  the  calculated  value  of  6*  by 
6*,  and  the  desired  value  by  6*,  we  form 


0[Pz(x)]  = 6*  - c*  (3.1.13) 

We  obtain  P2+^(x)  from  the  expression  given  by  (3.1.9).  To  find  the  deriva- 
tive of  ^ with  respect  to  P^  from  (3.1.13),  we  first  write  (3.1.4)  as 


6* 


(3.1.14) 


10 


where 


oo 

= J (c  — Ll)dr 


(3 


Differentiating  (3.1.14)  with  respect  to  p2  and  using  (3.1.13),  we  get 


^2  /R^ 


3A  A 1 3ue 


3p2  7 ^n  3p2 

e 


(3 


From  (3.1.5)  and  (3.1.15)  it  follows  that 


3U 


e _ n-1 


4a 


n 


3P2  Ue  ^pn+(l/?)  _ 


iff  = -F<0 


(3 


where 


F = 


_ 3f 


aPT 


f 

l 


.1.15) 


.1.16) 


.1.17) 


.1.18) 
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IV.  SOLUTION  OF  THE  GOVERNING  EQUATIONS  FOR  THE  STANDARD  PROBLEM 


4.1  Numerical  Formulation 

We  use  a very  efficient  and  accurate  numerical  method  to  solve  the 
governing  equations.  This  is  a two-point  finite  difference  method  developed 
by  H.  B.  Kel ler^ and  appl ied  to  the  boundary-! ayer  equations  by  Keller  and 
Cebeci  (see,  for  example,  references  16,  17). 


According  to  this  method  we  introduce  new  independent  variables  u(x,n), 
v(x,n),  t(x,n)  so  that  (2.3.4)  and  (2.3.b)  can  be  written  as  a first-order 
system 


f'  = u 

(4.1, 

.la) 

U ' = V 

(4.1, 

.lb) 

g'  = t 

(4.1, 

• lc) 

(bv) ' 

+ P1  fv  + P2  (c  - 

Z).x 

(4.1 

• Id) 

(et 

+ duv) 1 + P-j  ft  = x 

(u|a 

\ 3X 

-‘If) 

(4.1 

■ le) 

On  the  net  rectangle  shown  in  Figure  1 we  denote  the  net  points  by 


*o 

- 0, 

xn  = Vl  * kn' 

'i=l,2 N . 

no 

= o, 

"0  = V«  * V 

j = 1,  2,  • • 0: 

(4.1.2) 

nj  = 

Here  the 

net 

spacings  kn  and  h^ 

are  completely  arbitrary. 

The  quantities 

(f,  u,  v,  g,  t)  at  points  (xn,  rf j ) of  the  net  are  approximated  by  net  func- 
tions denoted  by  (f",  u",  v1?,  g",  t'J) . We  also  employ  the  notation  for  points 

V V J J J 

and  quantities  midway  between  net  points  and  for  any  net  function  q.: 


V',  5 2 (xn  * Vl>- 


"jV?  (nj  * nj-l> 


xn  - 


n-1 


n 1 # n , n . 
qj-'s  = qj-i> 


(4.1.3) 
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Fig.  1.  Net  rectangle  for  difference  approximations 


The  difference  equations  that  are  to  approximate  (4.1.1)  are  formulated  by  con- 
sidering one  mesh  rectangle  as  in  Fig.  1.  We  approximate  (4.1.1a,b,c)  using 
centered  difference  quotients  and  average  about  the  midpoint  (xn,n.  ;)  of  the 
segment  PjP2  as  follows: 

’’j'  (fj  - fj-l>  ' u"-'5  <4-’-4a) 

<uj  • uj-l>  = vj-b  (4.1 .4b) 

hj'  ‘ 9j-l!  = (4.1.4c) 


Similarly  (4.1.1d,e)  are  approximated  by  centering  about  the  midpoint 
(xn-*s»nj-Js)  0f  the  rectan^1e  P1P2P3P4 


<""’3  - b3-1vj-l>  * + - 


<P2  * 


/ n-l,n  ^n-1  n > 
+ an  v j -1/ j -h  ~ 


r""1  - P"  cn  , 


J-‘; 


2 \W, 


(h  . 1 ,4d) 
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k <«3*S  - * (p" * - “n  Mk  * vl*k 


n-1,  n . ..n-Un  tn-l ,n 
yj-H  j-*s  j-’s  J-'s  J-H 


Tj'’1 


(4.1.4e) 


where 


*k  • »„bn  - <u‘\n-]- H’  »TK'  - + p!',(f< 


+ p 


n-l 


*33  - <“2> 


n-l  1 1 

j-’sjl 


(4. 1 .5a) 


Tn-1 

T - l = a 

j-^  n 


<«>33  - 


n-1 

+ p1  [rth-h 


- 2 


[<*•>■]  33- 


hj'  (ej-'t-1 


en_1  tn-k 

ej-l  li-V 


(4.1.5b) 


Equations  (4.1.4)  are  imposed  for  j = 1,2,...,J.  The  boundary 
conditions  (2.3.7)  and  (2.3.8)  yield,  at  x = xp , 

i-H  a»  n n n ^ n n . n _ + n n i i * -i  r\ 

fo  ■ °*  uo=0*  UJ  = 1*  90  = gw  or  lo  " V 9J  ' K (4-1>6) 


4.2  Solution  of  the  Difference  Equations 

If  we  assume  (f^_1,  u1?'^ , , t"-1)  to  be  known  for  0 < j <_  J , 

J J J J v) 

then  (4.1.4)  for  1 j <_  J and  the  boundary  conditions  (4.1.6)  yield  a 

nonlinear  algebraic  system  of  5J+5  equations  ir.  ..s  many  unknowns 

(f!j,  u"  v"  g1?,  t1?).  The  system  can  be  solved  very  effectively  by  using 
J J 0 J 0 (i)  (i)  (i)  (i)  (i) 

Newton's  method,  We  introduce  iterates  [f'  uk  , vv.  .9^-  » , 

J J J J w 

i = 0,1,2,...,  with  initial  values  for  the  specified  wall  temperature  g , say 

W 
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.0,  u;0>  = o.  a«0,'V  ^0)  = 

-r-  “f-r- 

= g<0)-».  ‘J0)- 


C- 


tj-1  (4.2.1) 


For  the  higher  order  iterates  we  set 

,('*')  . ,(0  * „<<♦'>  . „<0  * Ju<<>,  v<’4|>  . ,!’>  * i,<!> 

J J J J J J J J J 


<,<<*"  . gh)  t igh\  . t(.‘>  ♦ 4t(i) 

J J J 


Then  we  insert  these  expressions  in  place  of  [f.,  u . , vH,  ys,  t.]  in  (4.1.4) 
and  drop  the  terms  that  are  quadratic  in  iv^,  i5j^’ 


This  procedure  yields  the  following  linear  system: 


5fj  ' ifj-l  ’ 21  (Suj  + Suj-1>  * 

h, 

5uj  - Suj-i  • r (5,j  * svm>  • (r2>j-^ 

69j  • ‘9j-i  •Ji(6tjtuj-;)  * (r3>j-' 


kl>j  6vj  * (t2>j  S,j-1  * (C3>J  <fj  4 U4>j  <fj-l  4 <S>j  5"j 


4 (t6>j  4uj-l  = (r4>j-H 


(4.2.2a) 


(4.2.2b) 


(4.2.2c) 


(4. 2. 2d) 


(Bl)j  6tj  + <*2>j  6tj-l  + (63'>j  4fj  + <e4>j  5rj-l  + (t5>j  ‘9j  + (66}j  6yj-l 
+ (pj)j  su^  + 6uj-l  = ^r5^j-c,  (4.2.2e) 
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for  j = Here  for  convenience  we  have  dropped  the  superscripts 

i and  n.  The  coefficients  (tK)i  (k  = 1 to  6)  of  the  differenced  imentun 
^uation  are: 


<r-l  >j  4 4 r [(pi  * “„)fj  - “n'j-'l] 

<C*)J  ■ 'bj-l  * 21  [!P1  * 

U3>j  ' 21  [(P1  * “n)wj  4 Vj-'J 
('4>j  " / [(P1  4 “n,vj-l  4 Vj-t] 


(4.2.3 


('5>j  ‘ ‘hj(P2  4 “n,uj 


(f.6)j  = -hj(P2  ♦ an)Uj_, 


The  coefficients  (e  ) . (k  = 1 to  H)  of  the  differenced  energy  equation 

k J 
are 


h . 


(6l)j  = ej  + r [(P1  + an)fj  - anfT-l] 

(S2}j  = 'ej-l  + / [(P1  + an)fj-l  * anfj-i] 


h . 

( R ) . = — J- 
\~'i  J -i  o 


3'j  - Z-  [<PI  4 %Hj  4 %£[] 


Wd * [(Pi + 


(4.2.4 


(4.2.4) 


(e8}j  = • T an  (9j-l  • 9j-i) 


The  coefficients  (r^).  (k  = ! to  5)  crc: 

vJ 


iri>j  ■ 'YW  fj  * fo-i 
(r2>j  “ Vj-h  - uj  4 “j-l 

(r3>i  * Vj-H  - 9j  + S0-<  (4.2.5) 

(r4> j ■ hj  (Rj-i  - P2cj-H>-(bjvj  ' bj-lvj-l  * hJ  [(P1  + 

- (p2  + %><“VS  - “n  (-*;:!  + f":!  v"-\»]  1 

(,'5>o  • Vjli  ' [Vj  ■ + hj  !(P1  * °nl(n)J-'. 

r/  \ n— 1 , n-1  , ,n-l  .n-K  ill 

- «„  [(“Slj.),  - VAV  Sj-H  “j-1-.  V.  V.  ■ 'J 


The  boundary  conditions  for  the  specified  wall  temperature  become 

= °»  = 0»  6Un  = 0*  CU,  = °»  *q,  = 0 


(4.2.6) 


The  boundary  conditions  for  the  specified  heat  flux  become 

itQ  = 0,  6t0  = 0,  5uQ  = 0,  6Uj  = 0,  -pj  = 0 


(4.2.7) 


We  use  the  block  elimination  method  discussed  by  Isaacson  and  Keller'  ' 
to  solve  the  linear  system  (4.2.2),  and  (4.2.6)  or  (4.2.7)  depending  on  the  boundary 
conditions.  Fo>  completeness  we  present  the  block  elimination  method  for  our 
problem  in  Appendix  A. 
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V.  SOLUTION  OF  THE  GOVERNING  EQUATIONS  FOR  THE  INVERSE  PROBLEM 
5 . 1 Var i_ajt  i onal  F guati ons 

In  order  to  calculate  a^/aP^  in  (3.1.10)  or  in  (3.1.16)  it  is  necessary 
to  know  F^  for  the  specified  case  or  F(n  )[  af/aPgO.  )]  for  the 

specified  6*  case,  for  this  reason  v/e  take  the  derivative  of  (4.1.4)  with 
respect  to  p£.  This  leads  to  the  following  linear  difference  equations,  known 
as  the  variational  equations  for  (4,1.4): 


- f;  o 

= U • , 

(5. 

1.1a) 

j-  • ; 

J —2 

-•w 

■ *!-■. 

(5. 

1.1b) 

(5. 

1.1c) 

^ * 1 <F!  * «„>  [W  * 

* 5-  - («*)?..,  - (*  * *„)  (W  * 

* 7 (S.Od) 

h-'(enTn  - e"  ,Tn  ,)  ♦ 1 (P?  ♦ , )kT"  * tnFn  * f"  ,Tn  . + F"  ,1 
J ' J J J-l  J-l7  2 ' 1 n'L  j j j J 3-1  3-1  3-1  J-1J 

, 1 >n  an  f nrn  , n..n  , n „n  , n ..n  , n-l/rn  , _n  \ 

+ -k-  (ft)  . i — -s-iu-G.  + g.U.  + u . ,G.  , + a.  ,U.  , + u . .(G.  +•  G . 

2 'j-1^  2 t J J yJ  J J-l  J-l  yJ'l  J-l  J-V  J J-l' 


( 5 . 1 . 1 e ) 


Equations  (5.1.1a  to  c)  can  be  written  a: 

h. 


Fj  - F"-i  - r (uj  + "j-i*  5 0 

"on  - ujn.i  - / (hn  + vj"-v»  ■ 0 


Equations  (5.1. Id,  e)  can  be  rearranged  in  a form  similar  to  (4. 2. 2d,  e),  that 
is: 


(5.1. 2d) 


(^T  )jTj  + (02)Tj_!  + (63)jF-  + (64)jF5_1  + ♦ (B7)jUj 


+ - 0 


(5.1.2o) 


Mere 


and  the  coefficients  ^ and 
and  (4.2.4).  The  coefficient 


8^  are  the  same  as  those  given  by  (4.2.3) 

(ra).  ] is  defined  by 
J— 5 


+ 


Similarly  the  boundary  conditions  (4.1.6)  become 


(5.1.3) 


Fn  - 0,  Un  = 0,  U,  = 0,  g"  * 0 or  T?  = 0,  g"  = 0 

(5.1.4) 

As  pointed  out  in  Ref.  12,  an  alternative  sec  of  variational  equations 

can  be  obtained  by  first  taking  the  derivative  of  (4.1.1)  and  then  differencing 

the  resulting  equations.  However,  this  procedure  does  not  necessarily  yield 

a good  approximation  to  the  desired  derivative,  V = .'v/'Pj  or  F - jf/,^. 

In  the  limit,  as  h.  and  k ->-0,  both  procedures  yield  the  some  result. 

0 n 

But  as  is  reported  in  Ref.  12,  for  the  actual  numerical  calculations,  the 
present  procedure  gives  precisely  the  derivatives  required  for  Newton's  method 
while  the  other  procedure  may  not.  This  is,  in  fact,  one  of  the  basic  differ- 
ences between  what  is  sometimes  called  ' quasi-linearization"  and  our  exact 
application  of  Newton's  method.  Thus  the  "quasi-linearized"  iterations  may  not 
converge  quadratical ly  (as  was  found  to  be  tne  case  in  references  12  and  17, 
but  our  present  iterates  do  show  super!  inea'"  convergence. 


The  system  (5.1.2),  (5.1.4)  again  forms  a block  tridiagonal  system  (with 
5x5  blocks)  that  is  easily  solved  by  the  block  elimination  method  described 
in  Apperiix  A. 
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VI.  RESULTS  FOR  STANDARD  AND  INVERSE  PROBLEMS 


6 . 1 Grid  Across  the  Boundary  Layer 

While  the  numerical  scheme  employed  here  is  a general  one  in  that  any 

type  of  grid  can  be  used  in  the  n-direction  (also  in  the  x-direction) , we 

(141 

have  chosen  a grid  previously  used  by  the  author  and  his  associates'  ' This 
grid  has  the  property  that  the  ratio  of  lengths  of  any  two  adjacent  intervals 
is  a constant,  that  is, 

hj  = Khj_]  (6.1.1) 

The  distance  to  the  j-th  n-line  is  given  ny  the  following  formula: 


nj  = h1  j = 1,  2,  3,  ...»  J i:  > 1 (6.1.2) 

There  are  two  parameters:  h^ , the  length  of  the  first  bn-step,  and  K,  the 
ratio  of  two  successive  steps.  The  total  number  of  points  J are  calculated 
by  the  following  formula: 


ln[l  + (K  — 1 )n  /h.] 

oc  | 

In  K 


(6.1.3) 


In  our  calculations  we  select  the  parameters  h-j  and  K and  calculate 
the  n . Several  runs  with  a different  number  of  points  across  the  boundary 

CO 

layer  showed  that  (see  ref.  14  ) approximately  30  to  40  points  are  sufficient 

for  turbulent  flows.  A typical  value  of  h,  is  0.01,  provided  that  the 

1 7 

Reynolds  number  is  not  very  large,  say  less  than  10  . At  higher  Reynolds 
numbers,  it  is  better  to  use  a smaller  value  of  h-j , say  h^  = 0.005.  Figure 
2 shows  the  values  of  K for  various  values  of  (n  /h,)  x 10~^  and  J. 

OD  | 

From  this  figure  we  choose  the  value  of  K as  follows: 

Let  us  assume  that  we  want  to  take  a maximum  of  40  points  across  the 
boundary  layer.  If  the  Reynolds  number  is  not  very  large,  an  estimate  for  the 
maximum  value  of  = 50  is  sufficient.  Then  taking  h-j  - 0.01,  the  ratio 
of  (n  /h,)  x 10‘?  is  50.  Thus  from  figure  2,  J = 40,  K is  approximately 

1.19 


20 


Figure  2.  Variable-grid  parameter  for  qiven  step  size  and  boundary-laye 
thi ckness, 


6 . 2 Results  for  the  Standard  Problem 

To  test  the  method  for  the  standard  problem  we  made  calculations  for  the 

(19) 

experimental  data  of  Lewis,  et  al  ■ , which  consists  of  compressible  adiahar.-'c 

turbulent  boundary  layers  in  adverse  and  favorable  pressure  gradients.  The 
results  are  shown  in  Fig.  3.  The  calculations  were  started  by  matching  a zero- 
pressure-gradient  profile  (Rd  = 4870)  at  x = 11.5  in.  downstream  of  the  leading 
edge  of  the  model.  Then  the  experimental  Mach-number  distribution  shown  in  figure 
3a  was  used  to  compute  the  rest  of  the  flow.  In  general  the  calculated  velocity 
profiles,  local  skin-friction  and  momentum-thickness-Reynolds-number  values  are  in 
good  agreement  with  experiment.  We  should  point  out  here  that  the  experimental 
skin-friction  values  were  obtained  by  Stanton  tube  and  were  not  deduced  from  the 
experimental  velocity  profiles 
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Figure  3.  Comparison  of  calculated  and  experimental  results  for  the  experimental 
data  of  Lewis  et  al.  (a)  Velocity  profiles  and  external  Mach  number 
distribution,  (b)  Local  skin-friction  coefficient  Cf  and  Reynolds 
number  R«  distribution. 


6.3  Results  for  the  Inverse  Problem 

To  test  our  method  for  the  inverse  problem  for  the  case  of  specified  cf-, 
we  made  calculations  for  incompressible  turbulent  boundary  layers,  and  checked 
the  results  obtained  earlier  in  another  study.  We  have  chosen  two  experimental 
incompressible  flows  from  the  data  reported  in  the  Stanford  Conference  on 
Computation  of  Turbulent  Boundary  Layers The  flows  we  have  considered 
are  known  as  1300,  5200  and  5300  in  that  conference. 

The  flow  1300  corresponds  to  an  accelerating  flow.  The  experimental  data 
is  due  to  Ludwieg  and  Tillmann.  Flows  5200  and  5300  correspond  to  decelerating 
flows  measured  by  Stratford.  They  differ  from  those  more  common  decelerating 
turbulent  flows  in  that  they  have  a negligible  skin  friction.  Thus,  they  are 
on  the  verge  of  separating.  For  this  reason  it  is  a very  severe  test  for  a 
numerical  method  and  for  exploring  the  accuracy  of  the  eddy-viscosity  formulas. 

In  the  1968  Stanford  Conference,  of  the  investigators  who  used  differential 
methods,  only  one  computed  5200  and  none  has  computed  5300.  The  accuracy  of 
computing  these  flows  is  also  important  in  many  design  problems  as  was  discussed 
in  the  Introduction.  The  design  of  the  Liebeck  airfoils  discussed  in  refs.  5 
and  6,  for  example,  is  based  on  the  results  of  boundary  layer  calculations  for 
a flow  with  vanishing  skin  friction. 

!n  making  these  computations  we  have  first  considered  the  standard  problem. 
That  is,  for  the  given  experimental  velocity  distribution  and  for  the  given 
initial  velocity  profiles  at  x = x0,  we  have  computed  the  velocity  profiles  and 
the  local  skin-friction  coefficient  at  each  specified  x-location.  Then  we  made 
the  calculations  for  the  inverse  problem.  We  specified  the  computed  local  skin- 
friction  coefficient  as  an  additional  boundary  condition  at  each  x-station  and 
computed  the  velocity  distribution  by  the  inverse  method.  We  have  thus  used  the 
computed  skin-friction  values,  rather  than  the  experimental  values,  as  a boundary 
condition.  Such  a procedure  is  necessary  because  a slight  error  in  the  experi- 
mental skin-friction  coefficient  will  severely  affect  the  computed  velocity 
distribution.  To  discuss  this  point  further,  let  us  consider  the  data  of 
Stratford,  either  5200  or  5300.  It  may  be  seen  from  the  skin-friction  p’ots 
in  figures  So  and  6b  that  the  experimental  values  of  cf  show  scatter;  in  an 
adverse  pressure  gradient  flow,  c^  should  either  stay  nearly  constant  or 
decrease.  If  one  uses  the  scattered  values  as  boundary  condition  and  computes 
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the  velocity  distribution,  one  would  get  slight  increases  and  decreases  in 
the  velocity  distribution  with  increasing  and  decreasing  cf)  respectively. 

The  computed  results  in  fig.  4 for  the  accelerating  flow  1300  show  very 
good  agreement  with  experimental  data.  This  indicates  that  our  eddy  viscosity 
formulation  is  quite  satisfactory  for  this  flow.  On  the  other  hand,  the 
computed  results  in  figures  5 and  6 for  the  two  decelerating,  on-the-verge- 
of-separating  flows,  5200  and  5300,  are  not  satisfactory  at  all  although  the 
computed  results  at  the  beginning  of  these  two  flows  agree  well  with  experi- 
ment (see  the  velocity  profiles  at  x = 2.9075  for  both  5200  and  5300).  This 
is  probably  due  to  the  effect  of  the  strong  pressure  gradient  suddenly  imposed 
at  x 3.0  ft.  However,  the  present  method  does  not  break  down  as  almost 
all  other  numerical  methods  do  and  performs  extremely  well  even  under  such 
strong  pressure  gradients  at  5200  and  5300  contain. 
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Figure  6.  Comparison  of  calculated  and  experimental  results  for  the  flow  5300. 

(a)  Velocity  profiles  and  external  velocity  distribution.  * 
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Figure  6.  (b)  Shape  factor  H,  local  skin  friction  coefficient  c, 

Reynolds  number,  R„  distribution.  f 
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Appendix  A 

BLOCK  ELIMINATION  METHOD 


The  linear  system  (4.2.2),  (4.2.6)  or  (4.2.7)  depending  on  the  boundary 

conditions  can  be  solved  in  an  extremely  efficient  manner  as  described  by 

fl  8 1 

Isaacson  and  Keller^  ' since  it  has  a b^ock  tridiagonal  structure.  This  is 
not  obvious  and  to  clarify  the  solution  procedure  we  write  our  system  in 
matrix-vector  form.  There  are  many  ways  in  which  this  can  be  done.  They  are 
all  equivalent  and  merely  amount  to  different  permutations  of  the  equations  or 
of  the  unknowns  or  both.  Further,  the  boundary  conditions  (4.2.6)  or  (4.2.7) 
could  be  employed  to  eliminate  the  five  unknowns  and  thus  slightly  reduce  the 
order  of  the  system.  When  the  latter  is  done,  the  system  (4.2.2),  and  (4.2.6) 
for  the  case  of  specified  wail  temperature  can  be  written  as 


where 


and 


ri  are  vectors  denoted  by 
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Here  the  coefficient  matrix  A is 


A1  C1 


B2  A2  C2 


(A. 3) 
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and  A ■ , C,,  B.  are  5x5  matrices  given  by  (a,  = h./^) 
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The  system  (A.l)  can  be  solved  by  the  block  tridiagonal  factorization 
procedure  described  by  Isaacson  and  Keller.  According  to  this  procedure,  we 
first  seek  a factorization  of  the  form 


t\  = XM 


(A. 5) 


where  (with  I.  denoting  the  identity  matrices) 
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From  (A. 3)  and  (A. 5)  it  follows  that 


“i  5 V 

VV  Vm 
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Substituting  (A. 5)  into  (A.l)  we  get 
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Then  from  (A. 7a)  we  find  that  for  j = 1 

i ^ 

'11  = (C2)1  -Tc7T[  1 


- + (^5^1  + a!  ^3)1 


*21  = TeiT 


-Ce'iTT  [(87)l  * a’(B3)'] 


y31  = _al 


Ml 


Ml 


(A. 12) 


35 
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y52  = ^ y22 

From  (A. 7b),  we  can  show  that  for  j >_  2,  the  elements  are  a^  are: 
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(a3i?)j  = _1  + aj(Y52)j-1 

( a42 ) j = ^ 8^j  “ ^ 4^  j ( y32  > j-1  " (^2VY42)j-1 

*“52Jj  = “ (B4)j(Y32)j-l  ~ (B2VY52)j-l 

The  unknowns  y^  in  (A.  11)  for  i = 1,  2,  3,  4,  5 and 
determined  from  (A.7c).  For  k = 1,  y^  are: 
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For  k = 2,  y^2  are: 
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<v32Jj  = ■(all)j(Y12)j  “ (a12)j(Y22)j 
( y42 ) j = ^a2lVYl2^j  + (a22VY22V 
(Y52)j  = £(a31  VY12}j  + (a32)j'Y22)j  “ ^ aj 


(A. 15b) 


where  (ai).,  (b,).,  (a2).j,  (62 ) , and  a.  are  the  same  as  In  (A. 15a);  (c2)i 
is  different  and  is  given  by 

<c^)j  * («5>j  + (6l>jaj' 

To  determine  w.  we  use  (A. 10)  and  the  definition  of  in  (fj).  If  we 
denote  by 
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Thus  for  j = 1,  using  th*  definitions  of  , w-j  and  r^  from  (A.  16a) 


can  write 


38 


. in  11  li'l  1 1 1 ililL  1 III  n ll^i  iitl  ^llu.  .k  ;llM  i I JuLii  uol  LoUii  u J h n n ' LQ  .t:.  1 1 Ll.tin.  J i I M JkixiiHUHtJUiiu.  Ji^iulllJ  1 iii:u  Lf  neU  Jill  luIiUUiij^iililliufi 


Thus  for  j - 1,  using  the  definitions  of  , w1  and  r-j  from  (A. 16a)  we 
can  write 
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where  (a,),,  ( ^ 1 ) - , (a?).,  (b?),f  A • are  given  by  the  definitions  in  (A. 15a) 
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Finally  the  perturbation  quantities  denoted  by  t are  calculated  from 
and  (A. 9).  It  follows  that 
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For  j = 1 


6v0  = ^1^1  “ Y,lVul  “ ^12^15ql 
*0  = ^w2^1  ” *Y21  VU1  “ *Y22V91 

6fl  = (w3}l  “ (y31)16U1  “ (y32}1691  (a.23) 

<Sv^  = (w4)1  - (y4])^u]  ~ ^42^1691 

6tl  = (w5^1  “ ^ y51 ^ 1 6U1  “ ^y52^1691 


Equations  starting  from  (A. 2)  to  (A.23)  are  for  the  case  of  specified  wall 
temperature.  When  the  wall  heat  flux  is  specified  some  of  those  equations  in 
(A. 2)  to  (A.23)  change.  For  example  6q  and  A-j  in  (A. 2)  and  (A. 4)  become 
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These  new  definitions  lead  to  the  new  equations  for  for  j = 1,  and  for 

(w^)-]  with  i = 1,  2,  ...»  5 and  k = 1,  2.  The  equations  for  y • k are: 

< Y31 ) -j  = "al 


( Y21 ^1 


(vn) 


(b6}1  “ (6lVl 


iri  ' - (c1)1 

*Y41*1  = “*Y11*1  - al 


rr  - (b3)1(y31)1] 

(Vi 

(Cc)i  ~ ( C-3  ) 1 ( Y-J1  )i  + “ ' C (Cr)-i(yoi) 


5'1  " ' ^3'1 ' y31 '1  T “a 


8yr  y21  'l 


( Y51  }1  = "(y21  )lal 


1 


(y32),  = 0 
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( Yp2 ) 1 


" (B«),  -(B,),.;'  C(B5)'  * (e’)|a'1] 


(y52^1 


= + ^y22^1  )3«71 


(>42^ 


(r-p1  ~ (C-|  ^C7^1  (r-8)i^22^1-' 


^12*1  = ”^42^1 


The  equations  for  (w^  are: 


(w3 ) 7 ' (rl^ 


(w2)1 


TT  C(r5),  -(63)1(w3)1  + (B1)1(r3)1a-1] 


(b6^  ^Bl^lal 

^1^1  " 1^2 rj"-  (?1)1  ^rVl  “ ^3Mw3^1  + ^1  Mr2^1al^  “ 

(w4)1  = -(w1)1  - (r2)1a“1  (A*26) 

“ ■[ (w2^i  + ^r3 -la] 

The  rest  are  the  same  except  for  the  second  equation  in  (A.?3).  That  equation 
is  replaced  by 


% = ^1  (y21>16u1  “ ^y22^1691 


(A. 27) 
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